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Abstract. Hedging methods to mitigate the exposure of variable an- 
nuity products to market risks require the calculation of market risk sen- 
sitivities (or "Greeks"). The complex, path-dependent nature of these 
products means these sensitivities typically must be estimated by Monte 
Carlo simulation. Standard market practice is to measure such sensitiv- 
ities using a "bump and revalue" method. As well as requiring multiple 
valuations, such approaches can be unreliable for higher order Greeks, 
e.g., gamma. In this article we investigate alternative estimators imple- 
mented within an advanced economic scenario generator model, incor- 
porating stochastic interest-rates and stochastic equity volatility. The 
estimators can also be easily generalized to work with the addition of 
equity jumps in this model. 



1. Introduction 

Many issuers of variable annuity (VA) contracts employ a hedging strategy 
to mitigate some of the risk borne out of writing these products. The chal- 
lenge of estimating the sensitivities of VA liabilities to their key risk drivers 
is, therefore, an important issue to many practitioners. For the purposes 
of this article a VA contract is "any unit-linked or managed fund vehicle 
-which offers optional guarantee benefits as a choice for the customer" 
The complex nature of the guaranteed benefits included in these products 
means that numerical techniques must be used to calculate such sensitivities. 
Monte Carlo simulation is perhaps the most fiexible of these approaches. 

In this article we extend the standard approaches to estimating option 
price sensitivities by Monte Carlo simulation to the problem of estimating 
the sensitivities of the liability of a stylized VA product. Furthermore, we 
adapt these approaches to be compatible with a more realistic economic 
model than the Black-Scholes framework, incorporating stochastic volatility 
and interest-rates. The main contributions of this article are as follows. The 
likelihood ratio method, a standard approach for evaluating option price sen- 
sitivities by Monte Carlo simulation, is extended to work under our more re- 
alistic market model with stochastic volatility and stochastic interest-rates. 
This builds on work by Broadie and Kaya [H], who only consider stochastic 
volatility. We also provide a full derivation of the equations for the likeli- 
hood ratio weights, which cannot otherwise be found in the literature. Also, 



Date: October 21, 2011. 



the pathwise method is developed in the context of VA habihties, follow- 
ing a similar approach to Hobbs et. al. [10], but considering a more complex 
VA product and the more realistic economic model. Finally, the pathwise 
and the likelihood ratio approaches are then combined together to construct 
a new and efficient mixed estimator for the second-order gamma sensitiv- 
ity of the VA liability. This will be important to practitioners looking to 
hedge their VA exposures, as obtaining an accurate and unbiased gamma 
estimate will allow the insurer to adapt their hedging strategy portfolio to 
take account of the convexity of their liabilities with respect to key market 
exposures. 

The paper is organized as follows; In Section [2] a overview of the standard 
methods for estimating option price sensitivities by Monte Carlo simulation 
is given. Three standard approaches are reviewed: the bump and revalue 
method, which is the natural finite difference approach often used in practice; 
the pathwise method, which was developed in the context of option pricing 
by Brodie and Glasserman j2]; the likelihood ratio method, developed in the 
context of option pricing by Broadie and Glasserman [2] and Glasserman and 
Zhao [6]. Mixed hybrid estimators, introduced by Broadie and Glasserman 
[2j, which combine the latter two of these standard approaches together to 
construct an efficient estimator for second-order sensitivities, will also be 
reviewed. 

In Section [3] the likelihood ratio method is extended to the setting of a 
more sophisticated economic model. This is an original contribution and 
builds on the work of Broadie and Kaya [3]. In Section |4] the standard 
approaches of Section [2] and the extension of the likelihood ratio method 
in Section |3] are developed for VA liability sensitivities. The pathwise ap- 
proach for VA liabilities follows a similar approach to the article of Hobbs 
et. al. [10], except that the present article considers a more complex product 
and the stochastic model for volatility and interest-rates. The likelihood ra- 
tio method is then extended to the liability of our stylized VA product. Next, 
the mixed estimators are constructed for the VA liability gamma sensitivity. 

To conclude the article. Section [5] compares all the estimators developed 
for the stylized VA product in Section [4] in terms of numerical efficiency. 
The mixed gamma sensitivity estimator is found to be particularly efficient, 
which is appealing as this is the sensitivity for which the standard approach 
performs the poorest. 



2. Option Price Sensitivity Estimators 

The liability from a VA contract takes a form very similar to the payoff 
of an option (or the sum of the payoffs of a series of options). Thus, in 
attempting to develop efficient methods for determining these VA liability 
sensitivities it would seem sensible to look for guidance from Monte Carlo 
methods for estimating option price sensitivities which have been developed 
by other researchers in recent years. We now summarize three classes of 
approach for calculating the sensitivities of options by simulation. Firstly, 
we introduce a simple method commonly adopted by insurers to estimate 
VA liability Greeks; the "bump &: revalue" approach. 
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2.1. "Bump & Revalue" Approach. This term refers to the concept of 
simulating the cost of the option under some base scenarios for key risk- 
drivers and then again under 'bumped scenarios', i.e., with the sensitivity 
parameter increased by some small perturbation, say /S.9. This sensitivity 
parameter could be the initial equity index level for a VA product, say. 
This is just a forward difference estimate in the sensitivity parameter. If 
the function Y[6) gives the discounted payoff of an option, then the price 
of this option is then given by a{9) = E[y(0)] with respect to a pricing 
measure. Now let Yi{9) , . . . ,Yn{0) represent the discounted option payoff 
along simulation paths 1, . . . ,n, the estimator of the first-order sensitivity 
would then be given by 

^ji_ Y{e + M)-Y{e) 

AO 

for the chosen 'bump size' A9, where Y{9) is the average of Yi{6), . . . ,Yn{9). 
The expectation of this estimator is given by 

^ a(9 + A9)-aiO) 
L J A9 ' 

The variance will be significantly reduced if the same random number 
stream is used in calculating the 'bumped' and 'base' option prices, as oppose 
to using independent random number streams. The 'bump and revalue' 
approach can incur fairly large sampling errors, particularly with second- 
order sensitivity estimates. Therefore, alternative approaches for estimating 
sensitivities are required by practitioners. Glasserman [5] introduces two 
more sophisticated general approaches to determining the sensitivities of 
option prices. 

2.2. Pathwise Estimator. The first of these approaches is the pathwise 
derivative method and the estimator of the first-order sensitivity with re- 
spect to 9 is given as follows. We can find the derivative of a{9) = K[Y{9)] 
analytically along each simulation path using 

= ,i,„ y(o+''i-ym, p.^ 

h^o h 

If the interchanging of differentiation and taking expectations is justified, 
that is if 



E 



^Ely(«)l 



then ^ Y17=i ^/(^) is an unbiased estimator of a' (9). 

Example 1. As an example of this method, consider the delta of a call 
option under the Black-Scholes model (i.e., the sensitivity of the option price 
with respect to 5*0 the initial underlying asset value). This can, of course, be 
found analytically, but helps illustrate the method. The discounted payoff 
of a call option is given by 

Y = e"'"'^ max(5T - K, 0) 

where r, K, T and St are the risk-free rate, strike price, time to maturity and 
equity price at maturity, respectively. Under Geometric Brownian Motion 

S'r = S'oe(^-5-')^+-v^^ (2.2) 
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where Z is a standard normal variate. Applying the chain rule to differen- 
tiate Y with respect to 5*0 with all other parameters held fixed gives 



dSo dSr dSo 

With the term the derivative fails to exist at the strike price, however 
the event {St = K} has measure zero and hence = e~^^I{ST > K) 
almost surely, where I{A) represents the indicator function of event A. For 



the second-term. Equation 2.2 gives = Thus, the pathwise estimator 



for the Black-Scholes delta is 

— -e I{St>K)-. 

If we wish to find the Black-Scholes gamma, i.e., sensitivity of the delta 
to the initial asset price, we have to differentiate W = e~^'^ I{St > K) with 
respect to So, which is equivalent to estimating the delta of a digital option. 
Here W is differentiable with respect to Sq with probability one, and takes 
the value zero. However, in this case 



= E 



dW 
dSo 



and the pathwise estimator is inapplicable (and is generally inapplicable 
when the payoff is discontinuous or in estimating second-order derivatives). 
The change in E[VK] with a change in 5*0 is explained by the fact that it could 
cause St to become 'in-the-money'. Glasserman [5j gives some technical 
conditions for when the pathwise method can be applied, however a less 
rigorous 'rule of thumb' is that it can be used when the payoff is continuous 
in the parameter of interest. 

2.3. Likelihood Ratio Estimator. The second more sophisticated ap- 
proach to estimating option sensitivities introduced by Glasserman [5] is 
known as the likelihood ratio method (LRM). The LRM approach relies on 
differentiating probability densities rather than payoff functions. Thus it 
does not require smoothness in the payoff function, as was required in the 
pathwise derivative method. 

Suppose we have a discounted payoff Y expressed as a function f{X), 
where X is a m-dimensional vector of different asset prices (or alternatively, 
one asset price at multiple valuation dates). Then, assuming that X has 
a probability density g with parameter 9, taking the expected discounted 
payoff with respect to this density gives 

E[Y] = E[f{X{e))] = [ f{x)ge{x)dx. 

Now, similar to the pathwise derivative approach, we assume the order of 
differentiation and integration can be interchanged. Here, however, this is 
not such a strong assumption, as typically densities are smooth functions, 
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whereas payoff functions are not. This gives 

J^m ge{x) 

= E 



= E 



9e{X) 
f{X)^^Hgo{X)) 



Then, f{X)^ln{gg(X)) gives the likehhood ratio estimator for the sensi- 
tivity with respect to the parameter 0, and this estimator is unbiased. The 
term ^\n{gg{x)) is often known in statistics as the "score function" and 
in this context is often referred to as the hkehhood ratio weight, since it 
multiphes the discounted payoff function to give the sensitivity estimator. 

The hkehhood ratio method will still be applicable and robust in the case 
of options with discontinuous payoff functions (and estimating second-order 
sensitivities) as this approach looks to differentiate the probability density, 
rather than the payoff function. 

Example 2. To show how this approach works in practice let us again 
consider the problem of estimating the Black-Scholes delta. In this case the 
discounted payoff function is given by /{St) = e~^'^ max(S'( — 0) and the 
log-normal density function used to calculate St is given by 



9{x) 



1 jH^)-^r-\a^)T 



xg\IT V o\fT 

where ^(•) represents the standard normal density function. In this case the 
score function for the parameter is: 

d ln(f ) - (r- ia2)T 

Evaluating this function at St and multiplying by the discoTintcd payoff of 
the option gives the unbiased estimator of the Black-Scholes delta as: 

e-^^ meociST - K, 0) ' . 

As we simulate for St using the relationship 

this likelihood ratio estimator can be written 

_rT — ^ 



e-''' max{ST - K,0) ■ 



SoaVr 

where the term Z / {Soas/T) is known as the LRM weight. 

One interesting point to note is the form of the payoff can be replaced by 
any other payoff function to give an estimator for the delta corresponding 
to the option with that payoff. For a digital option under the Black-Scholes 
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model we just multiply the payoff function for this option with the LRM 
weight which we have already calculated, giving the estimator for delta as 

Z 



e-^'l{ST>K) ^^^^^ 

This is a general feature of the likelihood ratio method which makes it 
appealing as a technique for estimating option sensitivities with different 
types of payoff in the same simulation run. Note that for reference, the 
LRM weight for the option gamma, found by an analogous process of dif- 
ferentiating the density geix) twice with respect to 9, can be easily shown 
to be z^-Z'^Vt-i 

For path-dependent options the LRM weights have a similar structure 
to the equivalent European option LRM weight for the required sensitiv- 
ity, except they will only use information along the path up until the first 
valuation date. The following Example demonstrates this. 

Example 3. The analysis thus far has studied how the LRM could be 
developed for European-type payoffs, however many VA products on offer in 
the markets will exhibit some degree of path dependency in their liabilities. 
With this in mind let us turn our attention to studying how the LRM extends 
to path-dependent options. 

Following Glasserman |5j, let us consider estimating the sensitivities of 
an arithmetic Asian call option with m valuation dates. The payoff on this 
option at maturity is then given by 

^ m 

Y = f{Si, . . . , Sm) = max{S — K,0), where S= — St- 

m ^-^ 

t=i 

Using the Markov property of Geometric Brownian motion, the underlying 
density for the equity asset path can be factorized as 

g{xi, . . .,Xm) = gi{xi\So)g2{x2\xi) ■ ■ ■ gmiXm\Xm-l) 

where gj{xj\xj-i) is the transition density from time to tj, i.e., 

9j{xj\xj-i) = ,/ =^{Cj{xj\xj-i)) 



(^j(Xj\Xj-l) 



log{xj/xj-i) - (r - a^/2){tj - tj.i) 



Suppose we wish to get an LRM estimator for the delta of the Asian 
option. From the above factorization it is clear that is a parameter of the 
first factor, gi{xi\So), only. This means we can express the score function, 
corresponding to the delta sensitivity, as: 

d\oggiS{ti),...,Sitm) ^ aioggi(5(ti)|So) ^ CiiSiti)\So) ^ Z, 

where Zi = Ci(S'(ti)|S'o) is the Gaussian increment which takes us from time 
zero to time ti. Likewise, the LRM estimator of the gamma sensitivity has 
a score function component which only relies on the equity asset path out 
to the first valuation date. 
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2.4. Mixed Estimators for Second-order Sensitivities. In estimating 
the second-order sensitivities of the option price with respect to some param- 
eter, for example the gamma sensitivity, one can also combine the pathwise 
and LRM approaches to create hybrid mixed estimators. 

Example 4. Again let us consider estimating the gamma Greek of a Eu- 
ropean call option under the Black-Scholes model. Firstly, by applying the 
pathwise approach to the LRM delta estimator, we obtain 

dSo\ SoaVrJ 
= e-'^I{ST>K)K ^ 



slaVr' 

Alternatively applying the LRM to the pathwise delta estimator gives the 
PW-LR mixed gamma estimator. The pathwise delta estimator has both 
functional and distributional dependence on Sq. To capture the distribu- 
tional dependence, the pathwise estimator is multiplied by the LRM weight. 
For the functional dependence another pathwise derivative is taken. This 
gives 



Sq Vcr\/T 



These mixed estimators typically give smaller standard errors for the gamma 
of European options, than the bump and revalue or pure LRM approaches 



3. Conditional Likelihood Ratio Estimator 

In order to use the likelihood ratio method, or the mixed hybrid second- 
order sensitivity estimators, one must be able to determine the probability 
density function of the underlying asset returns in order to derive the score 
function. With a Black-Scholes model, this score function is easily found in 
closed-form, but it is well known that this model does not give a realistic 
description of true market dynamics (see, e.g., [1]). There is, however, an el- 
egant approach which allows us to incorporate stochastic volatility /interest- 
rates and jumps in the equity returns, whilst still being able to utilize the 
tractability of the Black-Scholes model. 

The conditional likelihood ratio method (CLRM) was introduced by Broadie 
and Kaya jS] for the Heston stochastic volatility model [8j. In this article 
we extend this further to also incorporate stochastic interest-rates through 
the Cox-Ingersoll-Ross (CIR) model. Stochastic jumps can also be added to 
the equity process, but have been omitted here for ease of illustration of the 
method. See Broadie and Kaya j3] for the structure of the likelihood ratio 
weights when jumps are present, but stochastic interest-rates are not. In 
full, the system of Stochastic Differential Equations which govern our asset 
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price dynamics is: 

dSt = nStdt + ^/VtStdWf 

dVt = Kv{Ov - Vt)dt + avV^tdW^ (3.1) 

drt = Kr{9r - rt)dt + ar^/ndWl . 

Thus, the variance and interest rate both follow the same stochastic process 
(but calibrated with different parameters). The Brownian motions for the 
asset, volatility and interest rate processes are correlated as follows: 

corr(W-/, WY) = ps,v, corr(Wf , W^) = ps,r and corr(W,^, W^) = pv,r. 

In order to construct this correlation structure among the normally dis- 
tributed risk-drivers we employ a Cholesky decomposition. Let (^1,^2,-^3) 
be independent standard normal variates. Then standardized increments 
{Zv, Zr, Zs) for the variance, interest-rate and equity processes are con- 
structed by setting 

{Z^,Zr,Zsy = A{Zi,Z2,Z3). (3.2) 

Here A is the lower-triangular matrix satisfying AA' = p, where p is the 
correlation matrix constructed from ps,v, ps,r and pv,r- The correlation 
matrix p is assumed to be positive-definite. For completeness, the matrix A 
is given by: 



A 



f 1 \ 

PV,r \j^-p\r 

PS,r-pS,VPV,T L 9 {pS,r-PS,VPV,r)'^ 



The random variable Zr will now have correlation py^r with Zy and the 
random normal variate Zs will have correlation psy with Zy and corre- 
lation ps^r with Zr- This completes our Cholesky decomposition for the 
correlation structure of the three risk-drivers out to the first time-step. Now 
by conditioning on a realization of the variance and interest-rate processes, 
the asset returns become log-normal. One can then use the same form of 
LRM weights as for the Black-Scholes models, but with the addition of a 
couple of extra factors to account for the conditional information. Let us 
explore this idea further. Using the Cholesky decomposition the asset price 
dynamics can now be expressed as: 

^ = ndt + VVtdwf 

= ndt + vA^l^asidWfi + aaadW/^ + mzdW}^ 

where dWl^ , dW^^ and dW^'^ are three independent Brownian motions, 
corresponding to Zi, Z2 and Z3 in our discretization and aij is the element 
in row i and column j of matrix A. This construction for dSt expresses 
the asset price dynamics in terms of all the risk-drivers in the system of 
SDEs, which together dictate the behaviour of the asset returns. In theory 
an interest rate process with a larger number of risk-factors could be used. 
One would employ a Cholesky decomposition to correlate all the risk-drivers 



as required and express the asset price dynamics, dSt, in terms of all these 
other risk-drivers. One would however do this using a numerical Cholesky 
decomposition program to find the coefficients. 

To proceed towards a conditional log-normal representation of the asset 
returns, we need an expression for the stochastic process representing the 
logarithm of the returns, Xt- This expression for dXt, derived in Appendix 
A using Ito's formula, is given as 



with 



dXt = ndt + dYt - ^al^dt + V^ass dW^^' 



dYt = -^Vtal.dt - ^Vtal^dt + v^asidiy/^ + y^asadW^/^ 



Using this result for dXt, an expression for St in terms of the initial asset 
price 5*0 can now be found as: 

St = So eMYr) exp (^^ £ ndt - ^aj,^ £ Vtdt + aas^ £ v^dW/^) 
where Yt is given by 

T 7^ 



2 

By defining 



er = exp (Ft), (3.3) 
= rVtdt, (3.4) 




TT = ^ I ndt, (3.5) 



T 



'0 



the price of a European call option under the Heston SV with CIR interest 
rate model is given by C^^ {Sq^t, K,aT,rT,T), where C^^ is the Black- 
Scholes analytical formula for the price of a European call option. Also the 
CLRM weights can be determined using these adjusted expressions for Sq, 



a and r with the Black-Scholes score function, for example Equation 2.4 
for the delta sensitivity. This then gives us a means to obtain Monte Carlo 
simulation based estimators of the sensitivities of options under the proposed 
stochastic models for the underlying equity asset and risk-free interest rate. 
Naturally, the simulation must be performed in a conditional structure. 



i.e., the expectation is taken as E[/3] = E 
set-up is illustrated in Figure [T| 



E 



Vi,...,Vn,n, 



This 



4. Variable Annuity Liability Sensitivities 

4.1. Example Variable Annuity Product. The methods for estimating 
option price sensitivities will now be applied to the problem of estimating 
the sensitivity of the liabilities on a stylized variable annuity product. The 
idea behind this stylized example product is that it should be simple enough 
both in terms of tractability and ease of exposition, yet retain some of the 
key features which make these products so popular in many markets. It 
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should also have liabilities which are path-dependent, as is typical of many 
VA products on the market. This example product, based on a product from 
[TT] , will be of the Guaranteed Minimum Withdrawal Benefit (GMWB) type, 
where the policyholder is entitled to annual withdrawals from the underlying 
fund throughout the product lifetime, even if poor returns means that this 
fund diminishes to zero (at which point the insurance company must provide 
this income out its own reserves). 

More detail of this VA product and the policyholder will now be given. 
Firstly the contract owner is a 65 year old male and from the start of the 
contract has the 'guarantee rider' activated, i.e., this GMWB option initially 
sits on the pensions contract from the policyholder's 65th birthday. He 
will pay an extra 1% of the guarantee base (which is introduced below) in 
charges as a result of this option being activated, on top of the annual fund 
management charge which is deducted from the VA fund. If the policyholder 
wishes to turn off this option, he would be entitled to do so and this would 
cancel this extra guarantee charge. We shall refer to this as a customer 
'lapse'. The modelling of the rate of policyholder lapse is by no means a 
trivial issue and is worthy of study in its own right. We will assume a 
constant level of policyholder lapse of 4%, however a more realistic dynamic 
model of lapsation, based on fund levels and market conditions, could be 
employed without. This is a topic for further research. 

Under the GMWB contract the policyholder is guaranteed to receive in- 
come at the level of w% of the Guarantee Base each year after his 65th 
birthday, until he turns off the guarantee 'rider' or dies. This Guarantee 
Base is initially set at the amount of the policyholder premium, but can 
increase in value with an increasing VA fund level for the first a years af- 
ter annuitization. After this window passes, the Guarantee Base will then 
remain at the same level for the remainder of the product's lifetime. This 
ratchet feature, where the Guarantee Base steps up to the Fund Value if 
this is greater at the rebalancing date, is capped at a maximum year-on- 
year increase in the Guarantee Base of 15%. This will become clearer when 
the cashflows are described mathematically shortly. 
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The underlying VA fund, which is initially funded by the policyholder 
premium, is invested wholly in a single equity index with dynamics which 
are governed via the Heston-CIR model, given in Equations |3.1[ Possible 
extensions of this analysis include investing in a mixture of equities and 
bonds, or in a portfolio of two different equity indices. 

To model the cashflows on this policy mathematically let us define the 
Fund Value, Guarantee Base and policyholder income level at year t after 
annuitization by Ft, Gt and It respectively. Also let Rt denote the return 
from the equity index from year t — 1 to year t minus the management fees, 
i.e., Rt = St/St-i — 77, where tj is the fund management charge, quoted as an 
annual percentage. Then we can track the level of Fund Value throughout 
the lifetime of the policy using the following equation: 

Ft = max {{Ft-i-It-i){l + Rt),0). (4.1) 

This expresses the Fund Value at year t in terms of the Fund Value and 
Income level at year t — 1. Naturally then, we must have starting Fund 
and Income levels to initiate this recursion. The initial Fund Value is just 
given by the policyholder premium at annuitization, P (in units of the initial 
equity level), multiplied by the initial equity index level: Fq = P ■ Sq. The 
Guarantee Base at the end of year t after annuitization can be expressed as 

Gt = lit < a) min max(Gt_i, F*), 1.15 x + I{t > a)Gt-i. (4.2) 

The income level the policyholder withdraws from the policy Fund Value at 
the end of year t is then given by It = {'w—n)Gt- Here, iy is a fixed parameter 
dictating the proportion of the Guarantee Base that is withdrawn by the 
policyholder at each annual re-balancing date and is the guarantee charge, 
taken as a percentage of the Guarantee Base each year. In our case we have 
w = 4% and = 1%. 

The Liability, or Guarantee Shortfall, the insurer faces from issuing this 
VA contract on the market, measured at annuitization, can be expressed as 

" T 

L = E ^Apr™max(lt-Ft,0) (4.3) 
. t=i 

where Dt is the factor to ensure the liability at each yearly withdrawal date 
is discounted back to annuitization, p'^^'^^ is the probability of the policy 
remaining in force until year t after annuitization (encompassing both the 
possibility of policyholder mortality and lapsation) and T is the maximum 
contract term. Clearly, the insurer only faces a liability under the GMWB 
contract when the policyholder income cannot be met by the VA fund level. 
This is captured by the max function in the above formula for L, which sets 
the summand to zero at rebalancing/ withdrawal dates where Ft > It- 

4.2. Pathwise VA Liability Estimator. We begin by developing a path- 
wise methodology for estimating the VA liability sensitivities, for the exam- 
ple VA product of the last section. This method, proposed for a simple VA 
product by Hobbs et.al. [10], is just the natural extension of the pathwise 
approach for option sensitivities to the case of a VA liability, L. The liabil- 
ity is analogous to a series of European options of increasing maturity, thus 
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the same limitations of this approach for European option price sensitivities 
will apply here. Assuming that interchanging the order of differentiation 
and taking the expectation is justified the derivative of the liability with 
respect to the initial equity asset price can be expressed as 



Apw = IE 



T r. 

^Dtpr^—m^x{It-Ft,0). 



(4.4) 



The derivative only acts on the third term of the summand, and we can 
express the derivative of this term as 

-„,„(7,-F,.0) = /(7.>F,).(^-^). (4.5) 

The problem of estimating the delta sensitivity of the VA liability is now 
one of estimating the derivative of the Fund Value, Ft, and income level, 
It, for each year, t, after annuitization. Appealing to the structure of the 
product's cashflows these derivatives must be calculated recursively. Using 
Equation |4.1| we can express the derivative of the year t Fund Value with 
respect to 5*0 as 

dFt ( fdFt^i dlt-i\, 



Similarly, using Equation |4.2| the derivative of the year t income level with 
respect to 5*0 for the example product can be expressed using 

^ = l{{Ft>Gt-i)r^{Ft<l.lbxGt-i)r^{t<a))"^ 

+l{{Ft > Gt-i) n {Ft > 1.15 X Gt-i) n{t< a)) x 1.15 x 

obo 

+/((F,<G._.)n(*<„))«|^ + /(*>a)^^ 



where An i? is the intersection of events A and B, and the derivative of the 

dGt 
dSo' 



income level is then given by = {w — fi) ■ 



These recursions progress forward annually through the lifetime of the 
policy, with the initial conditions at time zero (annuitization) of |^ = P, 

the policyholder premium, and = 0, since the first withdrawal occurs at 



the end of year one. Substituting these recursion values into Equation 4.4 
via Equation |4.5[ gives the pathwise delta estimator. 



4.3. Conditional LRM VA Liability Estimator. The second alterna- 
tive Monte Carlo simulation approach for estimating option Greeks proposed 
earlier in the report was the (conditional) likelihood ratio method (CLRM). 
We now outline how the CLRM is used for estimating the sensitivities of 
liabilities arising on the example VA product. 

Firstly, we condition on a realization of the Heston variance and CIR 
interest-rate processes, in the same way as was performed in Section [3| The 
only difference here is that we only condition out to the first valuation date 
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at the end of year one, since the VA contract UabiHties are path-dependent 
(recall Example 2 in Section 2.3). That is, we calculate 

where Yi is given by 

-Jail [' Vtdt - lal2 C Vtdt + asi C ^/VdW}' + 032 ^/VdWl^ 
^ Jo ^ Jo Jo Jo 

The integrals in the above terms can be approximated by simple numerical 

quadrature. Then the (implied) shock out to year one is then given by: 

^ log(gi/6>go) - (fi - gf /2) X 1 
Using this implied shock the delta sensitivity can then be determined using 
A^^^ = (^^r|^) xLiabihty = (^^r|^) x ^ Apf ™ max(/t-Fi, 0) 

and similarly for the gamma sensitivity. The square-root of one is shown 
in the formula to help make the approach clear; it comes from the fact the 
first cashflow occurs at the end of year one. Thus, this term need not equal 
one if the product paid out semi-annual withdrawals, say, in which case this 
term would be \/0.5. 

4.4. VA Liability Gamma Mixed Estimator. Finally, we apply the 
pathwise approach to the CLRM estimator for the delta of the VA product 
liability in order to derive a mixed gamma estimator, in a similar manner 
as can be done for a European option. This results in 

pLR-PW ^ 9 LRM 

dSo 

j;Apr™niax(/i-Fi,0)) 
t=i J 

■V r(j . p^ (dit dFt 

T 

^Apr™max(/t-Ft,0). 

All the terms in the above formula are already obtained in calculating the 
pathwise estimator of the liability delta (except the LRM weight which is 
easily obtained). 

5. Numerical Comparison of VA Liability Estimators 

Having derived the delta and gamma estimators for the example VA prod- 
uct, we now study how the accuracy of the approaches compare. Five test 
cases will be considered, labelled A-E in Table [T| which give different pa- 
rameter settings for the Heston and CIR processes and different correlations 
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between the normal shocks driving the variance, interest-rate and equity 
processes. In all these cases the contract has term T = 30 years and the 
ratchet term is the first 10 years of the product lifetime. The income drawn 
each year by the policyholder is 4% of the Guarantee Base and the initial 
policyholder premium is £10,000. For these tests two separate simulation 
set-ups were created; one for the bump and revalue and another for the 
pathwise and CLRM estimates. In the bump and revalue set-up 36,000 
simulation paths were used. The bump perturbation size was set at 0.5%. 
There is a trade-off to be made here as the smaller the perturbation chosen 
the less bias in the estimate, however reducing the perturbation size will 
increase the variance of the estimator - particularly for the gamma. In the 
CLRM/pathwise set-up 10,000 variance/interest-rate outer paths with 10 
equity paths for each of these outer realizations were used. In both frame- 
works 20 timesteps per year were used to discretize the Heston and CIR 
processes. By doing this the bump and revalue and CLRM/pathwise ap- 
proaches took approximately the same amount of computation time to run, 
giving a fair basis for comparison of the approaches. 



Case 








'^IR 


%R 




PS-V 


ps-m 


PY-IR 


A 


2 


0.04 


0.15 


0.4 


0.04 


0.1 


-0.7 


-0.3 


0.2 


B 


1 


0.04 


0.3 


0.4 


0.04 


0.1 


-0.7 


-0.3 


0.2 


C 


2 


0.04 


0.15 


0.2 


0.04 


0.2 


-0.7 


-0.3 


0.2 


D 


1 


0.04 


0.3 


0.2 


0.04 


0.2 


-0.7 


-0.3 


0.2 


E 


1 


0.04 


0.3 


0.2 


0.04 


0.2 


-0.9 


-0.3 


0.2 



Table 1. Model settings considered in tests. Heston SV 
parameters denoted by V subscript, CIR parameters by IR 
subscript. Vq = ^y, ro = ^jj^. 





Sim. Set-up 1 (B&R) 
Liab.(£) St.Err 


Sim. Set-up 2 (PW/LRM) 
Liab.(£) St.Err 


A 


105.57 


0.68 


104.65 


1.18 


B 


125.19 


0.98 


123.68 


1.78 


C 


157.40 


1.02 


155.39 


1.81 


D 


169.11 


1.27 


166.86 


2.31 


E 


175.39 


1.52 


172.26 


2.87 



Table 2. VA Liability Estimates: VA example product. 



In table [2] the estimates of the liability value under the bump and revalue 
and CLRM/pathwise approaches are given. This table shows that under 
both simulation set-ups the values of liability are generally consistent. The 
bump and revalue approach does give estimates with smaller standard errors 
though. This is due to the conditional simulation framework required by 
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the CLRM approach, which is less efficient than a standard simulation set- 
up in estimating the product liabilities. However, it is the estimation of 
the sensitivity of these liabilities to the key market risk drivers that is the 
difficult challenge facing insurers. 





B&R St. Err 


PW St. Err 


CLRM St. Err 


A 
B 
C 
D 
E 


-0.00763 0.00013 
-0.00390 0.00014 
-0.00812 0.00016 
-0.00384 0.00017 
-0.00215 0.00019 


-0.00734 0.00014 
-0.00362 0.00017 
-0.00766 0.00019 
-0.00351 0.00021 
-0.00176 0.00029 


-0.00781 0.00032 
-0.00402 0.00040 
-0.00814 0.00042 
-0.00378 0.00050 
-0.00211 0.00060 



Table 3. Delta estimates: VA example product. 





B&R 


St. Err 


CLRM 


St.Err 


PW-LR 


St.Err 


A 


4.12 


0.55 


3.89 


0.47 


4.42 


0.09 


B 


3.25 


0.53 


2.83 


0.89 


3.17 


0.10 


C 


4.25 


0.77 


4.52 


0.70 


5.27 


0.12 


D 


2.50 


0.75 


3.45 


1.21 


3.83 


0.13 


E 


3.13 


0.70 


1.52 


3.77 


2.67 


0.21 



Table 4. Gamma estimates: VA example product. All 



numbers above should be scaled by a factor of 1 x 10 



Therefore, let us look at how the approaches perform in estimating the 
delta sensitivity of the liability in Cases A-E. Table [3] and Figure [2] give 
the estimates for the delta sensitivities for each of the approaches discussed 
earlier. Note that for the bump and revalue estimate, a central difference, 
using the 'bumped up' and 'bumped down' simulation paths, rather than 
a forward difference is used. This can help minimize levels of bias in this 
estimator. These results show that the bump and revalue and pathwise ap- 
proaches give similar estimates and standard errors for the delta in each 
of the given Cases. This is expected, as the pathwise estimator is essen- 
tially the small perturbation limit of the bump and and revalue approach. 
The reason the estimators are not converging to a common value is that 
they are simulated under the different set-ups: the pathwise results were 
estimated under the conditional simulation set-up alongside the CLRM es- 
timator, whereas the bump and revalue is from the base run of the bump 
and revalue set-up. Clearly the CLRM delta estimator is not as efficient as 
either the bump and revalue or pathwise equivalent (though still giving val- 
ues consistent with these estimates) . This can be explained by the fact that 
this method does not make use of the payoff function, unlike the pathwise 
approach. This, however, means the LRM estimator can still be used for 
options where this is discontinuous (and for second-order sensitivities). The 
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Comparison of Delta Estimators 
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'1 
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♦ B'nR Delta Est's « Pathwise Delta Est's ■ CLRM Delta Est's 



Figure 2. 



Comparison of Gamma Estimators 
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Figure 3. 
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pathwise estimator appears to be as efficient as the bump and revalue, but 
of course does not require a perturbed simulation run and is unbiased. 

In Table [4] and Figure [3] the different estimators for the gamma sensitivity 
of the VA liability are compared. These results show that both the bump 
and revalue and CLRM approaches give estimates with fairly large standard 
errors for cases A-E, with the CLRM being particularly poor for Case E 
due to the method not being very robust to high levels of /Jg.y- The mixed 
PW-LR estimator gives estimates which are consistent with these two ap- 
proaches, however they yield a much smaller standard error for each of the 
Cases considered. 

Thus, by constructing a mixed estimator for the gamma sensitivity, we 
have managed to gain a significant amount of reduction in variance, com- 
pared to the bump and revalue approach. Furthermore, this method is 
unbiased, hence we do not encounter the problem in choosing the pertur- 
bation step-size, where there is a trade-off between trying to reduce bias, 
yet also keep the size of the standard error of the estimate small. This is 
consistent with the analysis of the mixed estimators applied to European 
options in Glasserman [5J. 

This significant improvement in efficiency in estimating second-order sen- 
sitivities could be of great benefit to practitioners in managing a hedging 
strategy for their VA books. With the simple bump and revalue approach it 
is difficult to get an accurate and unbiased estimate of the gamma sensitiv- 
ity. As a result insurers will find it difficult to appreciate how frequently to 
rebalance delta-hedge positions due to underlying market movements. With 
a much more accurate and unbiased gamma estimate, the insurer should 
have a much greater appreciation of the convexity in the liabilities in their 
VA books and can adapt their hedging strategy portfolio accordingly. 

6. Conclusion 

With the increasing popularity of VA products and the new Solvency II 
regulatory framework in Europe, employing an effective hedging strategy 
for mitigating the risks from marketing such products is a current challenge 
facing insurers. The recent financial crisis has demonstrated that under tur- 
bulent market conditions a hedging portfolio can require much more frequent 
rebalancing. The widely adopted bump and revalue approach for estimating 
the Greeks used in such a strategy has some shortcomings, particularly for 
second-order sensitivities, such as gamma. 

In this article some more advanced estimators for VA Greeks have been 
developed which are unbiased and do not require additional perturbed sim- 
ulation runs. The mixed estimator developed for the VA Liability gamma 
sensitivity has much smaller standard errors compared to the bump and 
revalue method. Furthermore, the bias-variance trade-off in choosing the 
perturbation size in the bump and revalue is avoided. The estimators are 
also comparatively more efficient as the required number of Greeks increases. 
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Appendix A. Derivation of Stochastic Process for Xf. 



In this Appendix, a derivation of the formula for the stochastic process 
representing the logarithm of the returns, Xt, is given. Recall, that by 
employing the Cholesky decomposition the asset price dynamics under the 
Heston-CIR model were expressed in Section [3] as 

^ = ndt + VVtdwf 

= ndt + v^l^agidWf ^ + a^2dWl' + a-^-^dW^ 

where dW^^ , dW^^ and dW^'^ are three independent Brownian motions. 

Using this formula the expression given in the article for the stochas- 
tic process representing the logarithm of the returns, Xt, can be found as 
follows. Let g{t,x) be such that Xt = g{t, St) = In St- Then, 

dg _ Q dg _ 1 d'^g _ 1 

dt ^ dx X dx'^ 



Using Ito's formula, dXt, the stochastic process for the logarithm of the 
asset returns under the Heston-CIR model, is then given by: 

dXt = ^{t,Xt)dt + ^{t,Xt)dSt + l^{t,Xt){dStf 
= ndt + ^/Vt{audWl^ + a32<iVF/' + aasdV^/^) 



1 1 



2 52 



(Stndt + St^/Vt(a^ldWl^ + a^2dWl^ + a^^dW^'^ 



ndt + ^A^a3l^iVF/l + v^a32diy/^ + ^V^assciVF/^ 

111 

--Vtalidt - -Vtal2dt - -Vtaj^dt 

ndt + dYt - ^al^dt + ^A^a33^i^V/^ 
with 

dYt = -^Vtal^dt - ^Vtal^dt + v^agidVF/^ + y^ta32dWt^\ 
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